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The interaction between two spherical colloidal particles with degenerate planar anchoring in a nematic media is studied by 
numerically minimizing the bulk Landau-de Gennes and surface energy using a finite element method. We find that the energy 
achieves its global minimum when the particles are in close contact and making an angle = 28° ± 2 with respect to the bulk 
nematic director, in agreement with the experiments. Although the quadrupolar structure of the director field is preserved in 
the majority of configurations, we show that for smaller orientation angles and at smaller inter-particle separations, the axial 
symmetry of the topological defect-pairs is continuously broken, resulting in the emergence of an attractive interaction. 
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1 Introduction 

Studying the behavior of colloidal particles in anisotropic flu- 
ids with long-range orientational ordering, such as nematic 
liquid crystals, has attracted a great attention in soft condensed 
matter physics^"—. The orientation order parameter of the fluid 
(e.g. director of nematic liquid crystal) is distorted from its 
uniform orientation in the bulk due to anchoring to the surface 
of the colloidal particles. These elastic distortions create topo- 
logical defects around the particles 7 and induce anisotropic 
long and short range interactions between the particles^. 

Depending on the colloidal material and its coating, the 
surrounding fluid may have a normal orientation (normal or 
homeotropic anchoring), or parallel orientation (planar an- 
choring) with respect to the colloidal surfaces. For normal 
anchoring, the orientation of the fluid is locally and uniquely 
determined on the colloidal surfaces. For planar anchoring, 
the orientation of the fluid is degenerate on the colloidal sur- 
faces and is determined by the global structure of the fluid. 
However, the director field in the surface of the particles is af- 
fected by environment for any finite anchoring, such a free- 
dom makes the theoretical investigations more complicated 
for planar anchoring. 

In case of a single colloidal particle, the particle-defect pair 
induces a dipolar or a quadrupolar long-range elastic distor- 
tion fiel d 10 ' 11 depending on the anchoring type. The long- 
range dipolar structure results from a satellite point defect, 
when the size of the particle is large compared to the co- 
herence length of the nematic fluid and the anchoring is nor- 
mal^. The quadrupolar configuration appears in both normal 
and planar anchorings. In the normal case, a disclination ring 



" Department of Physics, Sharif University of Technology, P. O. Box 14588- 
89694 Tehran, Iran. 

b Department of Physics, Harvard University, Cambridge, MA 02138, USA. 
c Nanosystem Research Institute, National Institute of Advanced Industrial 
Science and Technology (AIST), 1-1-1 Umezono, Tsukuha 305-8568, Japan. 
* ejtehadi@sharif.edu 



(saturn ring) encircles the particle, when the size of the parti- 
cle is small or the strength of anchoring is weak- 7 ^. In planar 
anchoring, the elastic distortions form two point defects (boo- 
jums) at the poles of the particles, aligned along the nematic 
direction- 7 ^. 

The more physically interesting configurations are achieved 
when there are many colloidal particles present in the medium. 
For large separations of particles, the defects around of each 
of the particles is independent of that of the other particles, 
and (anisotropic) interaction potential between them is deter- 
mined by the long-range orientational field of the fluid. In this 
regime, and in case of two particles separated by a distance 
d, the effective interaction potential between them is propor- 



tional to d~ or d for dipolar or quadrupolar defect con- 
figurations, respectively 1 i 3 i 7 i 8 . When the particles approach 
each other, the defect structures are distorted and the interac- 
tions deviate from the far-field dipolar-dipolar or quadrupolar- 
quadrupolar interactions—"—. 

Experimentally, the colloidal interactions in nematic liquid 
crystals are studied using optical— or magneto-optical tweez- 
er s 19 ' 20 . The medium-induced interactions play an essential 
role in the formation of chain— or crysta l 2 - 21 i 22 suspensions of 
the colloids (or droplets) in nematic liquid crystals. 

Smalyukh et ai— have measured the angular and the ra- 
dial components of the force between two particles with pla- 
nar anchoring in a nematic liquid crystal as a function of the 
inter-particle separation d, the angle between the bulk ne- 
matic director and the vector connecting the particles, (see 
Fig. [TJ. They observe deviation from the far-field theoreti- 
cal quadrupole-quadrupole interaction— when the objects are 
in the close-contact regime. They also find that the equilib- 
rium configuration corresponds to a very small separation of 
particles, d ~ 2R, and ~ 30°. 

In this paper, we study the fluid-induced interaction be- 
tween two spherical colloidal particles of radius R by numer- 
ically minimizing the sum of elastic Landau-de Gennes free 
energy of the bulk fluid— and the degenerate planar anchor- 
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ing surface energy introduced by Fournier et ai^. In par- 
ticular, we are interested in the regime of strong anchoring 
and large particles (R large compared to the coherence length 
of the fluid, £). Our main goal is to study the close-contact 
configurations for which no theoretical work has been done 
to our knowledge, though interesting physics is expected to 
emerge due to the strong interaction of the topological defects. 
We also aim to explain the experimentally observed angle of 
~ 30° at equilibrium. 

The paper is organized as follows. The theoretical model is 
explained in Section [2] We describe the details of our numer- 
ical approach in Section [3] We finally present and discuss the 
results in Section|4]and summarize our findings in Section|5] 



2 The Model 

The geometry of the studied system is schematically illus- 
trated in Fig.Q] We consider two identical spherical colloidal 
particles with radius R immersed in a 3D nematic cell. The ne- 
matic director is aligned along x-axis at the boundaries of the 
cell. We scale all the lengths with respect to the radius of the 
particles. The dimensions of the cell is L x = 15R, L y = \5R 
and L z = 6R. The centers of the particles are confined to the 
plane z = L z /2 and are separated by a center-to-center distance 
d. The line joining the center of the particles makes an angle 
8 with x-axis. The dimensions of the cell are chosen in a way 
to ensure that in all of the studied configurations, the distance 
between the boundaries and the particles is much larger than 
the coherence length and big enough so that the nematic di- 
rector distribution is not affected by the boundaries. The free 
energy of the system can be written as: 



F(d,Q) = F h {d,Q)+F s (d,Q)+F c _ c (d), 



(1) 



where F\,(d,Q) is bulk nematic fluid free energy , F s (d,Q) is 
the surface energy and F c - C (d) is the Van der Waals colloid- 
colloid interaction. The free energy functionals will be de- 
scribed in detail in the following sections. 

In realistic situations, the only appreciable effect of Van 
der Waals colloid-colloid interactions is to provide a short- 
range repulsion between the colloidal particles. Such effects 
become relevant only in the regime where the separation of the 
colloid surfaces approaches the atomic length-scales. In this 
study, we confine ourselves to the regime where the surface 
separations are larger than the nematic coherence length, i.e. 
d — 2R 3> % 1 A. We note that the surface separations can 
chosen to be appreciably smaller than the size of the particles 
in this regime. Therefore, due to the separation of scales in 
this regime, we ignore the Van der Waals interaction between 
the colloids in this study. 




Fig. 1 Schematic representation of the studied system. The nematic 
director is fixed in x direction at the boundaries and two particles are 
symmetrically placed in x-y plane with respect to the center of the 
box. 



2.1 The Nematic Order Parameter 

The nematic fluid is described by a local, 3x3, traceless 
and symmetric tensor order parameter, Qij = S(hihj — 8y/3), 
which can be specified by five independent components, 



Git Qn Q13 
Qn Q22 623 
Q 13 Q23 Q33 



(2) 



where 633 = —{Qn +Qi2)- The scalar order parameter, S, 
and the director orientation, h, are locally obtained by the 
largest eigenvalue of the tensor order parameter, A, max = IS, 
and its corresponding eigenvector, respectively. We note that 
working with a tensor order parameter and expanding the free 
energy in its terms allows the formation of biaxial order, which 
is a necessary ingredient for a realistic description of topolog- 
ical defects and their interactions. 



2.2 The Bulk Free Energy 

The bulk free energy of the nematic fluid is described well 
by the Landau-de Gennes model—, in which the free energy 
functional is expanded in powers of the tensor order parameter 
and its spatial derivatives: 

Fb = L dv {^ QuQji - f Q v Q J kQki + j (QijQjt) 2 



u 

2 



(3) 



where the indices refer to Cartesian coordinates, the summa- 
tion over repeated indices is assumed and Q. denotes the vol- 
ume occupied by the nematic liquid crystal. The first three 
terms are the Landau-de Gennes free energy which describe 
the bulk isotropic-nematic (IN) transition. The coefficients A, 
B, and C are the material-dependent parameters. 
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The derivative terms are the contribution of the elastic 
free energy in the nematic phase. The nematic elastic con- 
stants, L\ and L2, are related to the Frank elastic constants 
by Li = K twi;i( /2S 2 and L 2 = (^ sp i a y - K twist )/2S 2 = (K bend - 
^twist) /2S 2 . In this study we restrict ourselves to one-elastic 
constant approximation that means all the Frank elastic con- 
stants should be equivalent which leads to L2 = 0. 

To simplify calculations, we rescale the tensor order param- 
eter as q = (4B/3\/6Cy l Q, such that q t j = S{hihj - 8,7/3), 
where 5= (4B/3y/6C)~ 1 S. As a consequence, the dimension- 
less free energy becomes 

^ b = Jo& = J a \2 qijqji _ ~4~^J q j kqki + 4 WW 

(4) 

where f = C(45/3%/6C) 4 , X = 27 AC/85 2 is effective di- 
mensionless temperature, \ = ^J21L\C /85 2 is the nematic 
coherence length 25,26 . In these dimensionless units, the 
fluid undergoes a first-order isotropic -nematic transition at 
T = 1/8. The isotropic phase becomes unstable for x < 
. The scalar nematic order parameter in bulk is given by 

5 b = (3V6/16) (l + v/1 -64x/9) . Throughout this paper, 
the quantities appearing with a hat are rescaled with respect 
to the the radius of the particles, e.g. dV = /?~ 3 dV, ^ = R~ 1 ^ 
and d k = Rd k . Following the related previous studie s 15 i 25 i 26 ; 
the dimensionless temperature and length scales were set to 
x = (3^6 — 8)/12 and q = 0.03 respectively. This choice of 
parameters closely match the parameters of the widely used 
liquid crystal mesogen 5CB and results in formation of stable 
topological defect s 25 i 26 . 

2.3 The Surface Free Energy 

As mentioned in the introduction, the preferred anchoring can 
be modeled with much more eas e 27 ' 28 compared to the planar 
degenerate anchoring due to uniqueness of the orientation of 
nematic fluid in the vicinity of the colloidal surfaces. 

Fournier et al— have recently introduced a two-parameter 
surface energy functional that is bounded from below and as- 
sumes its minimum in the manifold of degenerate planar con- 
figurations. The surface energy consists of two terms, control- 
ling the planar anchoring and fixing the scalar order parameter 
on the surface. Since we expect the formation surface topo- 
logical (e.g. two surface defects of charge +1 jl in case of a 
single colloid), we relax the second constraint like previous 
studies 2 -, resulting in the following single parameter surface 
energy functional: 

F s = W [ dA (Qij - Qjj)(Qji - Qji), (5) 
Jdn 



Fig. 2 A typical tetrahedral mesh used in the FEM analysis. The 
mesh is finer near the colloidal surfaces to capture a more accurate 
representation of the curved surfaces and to provide increased 
numerical accuracy near the topological defects. 

where Qij = gy+55 i7 /3, Q± = {?>i k -ViV k )Q k ,(?>ij-ViVj) is 
the projection of Qij onto the tangent plane of the surface, and 
v is the normal to the surface. W is positive and controls the 
stiffness of anchoring. The surface energy can be written in 
dimensionless variables: 

~ Fa W f I I 

where w = 27WC/8B 2 , q tj = (45 fiVbC)- 1 Qij, qjj = (8 ik - 

ViVk)qici(&ij — V/V 7 -) and dA = R~ 2 dA. We chose w = 0.0156, 
which describes strong planar anchoring to the surface once 
we take the bulk free energy parameters into account. 

3 Numerical Minimization of the Free Energy 

We adopt a finite element method (FEM) approach to min- 
imize the free energy functional described in the preceding 
sections. The nematic cell was decomposed into tetrahedral 
elements by using the automatic mesh generator GmsliM. In 
order to capture a more accurate representation of the curved 
surfaces and to provide increased numerical accuracy near the 
topological defects which are expected to be formed on or in 
the vicinity of the surfaces, a finer mesh size of Lsms = 0.025/? 
was used near the shperical boundaries. We note that since the 
nematic coherence length is £, = 0.03/?, the physics of elastic 
deformations can be properly captured in the used mesh. The 
mesh size was increased to Llms = 0.25/? away from the parti- 
cles in order to reduce the computational cost (see Fig.|2]i. The 
tensor order parameter was linearly interpolated within each 
of elements in the evaluation of the free energy integrals. We 
note that linear interpolation is the simplest scheme that pre- 
serves the properties of Q as an order tensor. For each configu- 
ration of the particles, the total dimensionless free energy was 
minimized using a conjugate gradient (CG) method—, yield- 
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Fig. 3 The free energy landscape of the system of two particles as a 
function of the inter-particles distance, d/R, and 0, the angle 
between the line joining the center of particles and the far-field 
director. 

ing the Effective Potential Energy (EPE) U(d,Q). In order to 
accelerate the minimization procedure for any given (8, d), 
we used the relaxed director profile from the closest preceding 
configuration as the initial guess. The configuration space of 
the particles was scanned in the range 2.10/? < d < 3.5R and 
0° < 9 < 90° in steps of 0.05R and 2° respectively. For larger 
separations, 3.5R < d < 6.0R, we scan with larger 8 steps of 
10°. For each configuration, the minimization procedure was 
stopped when the relative free energy improvements dropped 
below 10~ 5 . We present and discuss the results in the next 
section. 

4 Results and Discussion 

Figure [3] shows the full equilibrium free energy landscape of 
the system. The smooth contour plot of the free energy land- 
scape is shown in FigHJfor better clarity. The normal to the 
contour lines specify the direction of the net force between the 
colloidal particles. The effective potential energy of the sys- 
tem of two colloids is also shown in Fig. [5] as a function of 
the particle separation d, for four different orientations. The 
figure shows that for = 90°, the particles repel each other in 
the whole range of inter-particle separations, while they attract 
each other at = 60° and = 30°, with a stronger attraction 
in the later case. This uniform attractive or repulsive behavior 
is destroyed for configurations with smaller angles. In partic- 
ular, in the case of = 0, the particles attract each other in the 
close-contact regime, d < 2.7 R, while they repel each other for 
the larger separations. We will show later that this peculiar be- 
havior is associated to the spontaneous broken axial symmetry 
of the defect pairs. 

In order to analyze the distance dependence of the ef- 
fective potential, we have fitted a two-parameter function 
cq + c\ (d/R) Cl to the plots of Fig. [5] The fits are shown in 




2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 
d/R 



Fig. 4 The contour plot of the EPE between two particles as a 
function of the inter-particles distance, d/R, and 6, the angle 
between the line joining the center of particles and the far-field 
director. 



Fig. |6]and correspond to exponents c-i = — 5.6 ± 0.2 and C2 = 
-5.7 ± 0.2 for = 0° and = 90° respectively. It is noticed 
that exponents are slightly larger in magnitude in comparison 
to the weak anchoring analytical analysis—, C2 6ak = —5. The 
exponent is c? = —2.5 ±0.4 for = 60°, which describes a 
deformation field of longer range in contrast to the weak an- 
choring theory. Finally, in the case = 30°, the exponent is in 
agreement with the analytical prediction, C2 = —5.0 ±0.2. 

Fig. [7] shows the effective potential energy as a function 
of angle for fixed particle-particle separations (d/R — 2.10, 
2.30, 2.50, 2.70, 2.90, 3.10, 3.30 and 3.50). It is noticed that 
each of the plots has a unique global minimum. Moreover, the 
free energy minimum, as well the orientation angle at which 
the minimum is achieved (0 m i n ), decreases monotonously as 
the particles approach each other. This behavior is shown 
clearly in the inset plot of Fig. [7] Extrapolating m ; n to the 
limit d ~ 2R, we find that the global minimum of the free en- 
ergy is achieved for = 28° ±2°. This result is consistent 
with the experimental results of Poulin et al.—, Smaylukh et 
al— and Kotar et al.—, who find the equilibrium angle to be 
~ 30°. 

The vector field of the net force between the particles can 
be calculated by taking the gradient of the effective potential 
energy. The net force field is shown in Fig. [8] and the configu- 
rations at which the radial (F r ) and tangential (Fq) components 
of the net force vanishes are indicated. It is easily noticed that 
the force field drives the system towards the configuration of 
minimum free energy, i.e. m ; n ~ 28° and d ~ 2R. The radial 
component of the net force is positive for > 60° and thus, the 
force is repulsive. The angle at which the radial component of 
net force changes sign depends on the inter-particle separa- 
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Fig. 5 Effective potential energy between two particles as a function 
of the inter-particles distance, d/R, for four different orientations 
6 = 0°, 30°, 60° and 90°. 
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Fig. 7 Effective potential energy as a function of angle 9 for 
different inter-particle separations d/R = 2.10, 2.30, 2.50, 2.70, 
2.90, 3.10, 3.30 and 3.50. The inset plot shows the angle at which 
the free energy achieves its minimum, 9 m ; n , as a function of 
inter-particle separation. The error bars indicate numerical 
uncertainty in the minimum position. The free energy assumes its 
global minimum at 6 m ; n = 28° ± 2° . 
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Fig. 6 Log-log plots of the effective potential energy vs. 
inter-particle separation shows an orientation-dependent asymptotic 
power-law behavior. Fitting a two-parameter function 
c + C] (d/Rf 2 to U(d,B), we get c 2 = -5.7 ±0.2, -5.0±0.2, 
-2.5±0.4 and -5.7±0.2 for 6 = 0, 30°, 60° and 90° respectively. 



tion and varies in the range 60° < 8 < 70° for the investigated 
configurations. Repulsive interaction is expected to show up 
when > 75° in quadrupolar approximation—. We associate 
this slight discrepancy to deviations from the quadrupolar ap- 
proximation. 

As it was mentioned earlier in this section, the inter-particle 
force exhibits a non-monotonous behavior as a function of 
inter-particle separation for small angles (0 < 15°). This be- 
havior can be seen in Figs.|5j|4]and[8] In particular, for = 0°, 
the net force is repulsive for large separations, while it be- 
comes attractive for ~ d/R < 2.7. To our knowledge, this pe- 
culiar behavior has not been studied theoretically elsewhere. 

Moreover, its experimental observation requires measure- 
ment the interaction force in separations smaller than those 
reported in the experimental paper by Smalyukh et a/.^i. 
We note that precise experimental measurement of this phe- 
nomenon can be carried out by fixing the orientation of the 
colloidal particles by using line optical tweezers^ 2 -. 

In order to gain insight into this phenomenon, we study 
the configuration of the topological defects at different inter- 
particle separations, as shown in Fig. [9] When the particles 
are well separated (d/R> 3), the boojum defects are aligned 
on the x-axis and the defect-particle pairs have a quadrupolar 
symmetry for each of the colloidal particles. The repulsion in 
this regime is thus associated to the head-to-head interaction 
of defects of equal charge +1/2 (Fig.|9^). When particles ap- 
proach each other, the continuous 0(2) symmetry of the de- 
fect pairs is continuously broken due to the strong repulsion 
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Fig. 8 The vector field of net force between two colloidal particles 
in the x-y plane. The center of one of the particles is fixed at the 
origin. F r and F$ denote the radial and tangential components of the 
force, respectively. 



between the approaching boojums, driving them away from 
the axis of symmetry (Fig. |9j)-c). The nematic director pro- 
file on the front and back of the colloidal particles is shown in 
Fig. [TOb andfTUb respectively, as seen along the x-axis. It is 
noticed that the approaching pair of boojums are driven away 
from the .Y-axis (Fig. [TOh), while the boojums on the back of 
the particles remain on the x-axis. The displacement of the ap- 
proaching pair of boojums induces an attraction between the 
colloidal particles due to the energetic tendency to reduce the 
volume of the distorted region between the particles. 

At larger angles (0 > 30°), the axial symmetry of the defect- 
pair on each of the particles is almost preserved in the whole 
range of inter-particle separations (see Fig. QT|and Fig. \l7\ . 
Therefore, the monotonous attractive or repulsive interaction 
at larger angles can be explained by merely taking into account 
the quadrupolar deformation field of each of the particles—. 

We finally remark that although the calculations were car- 
ried out for a specific choice of bulk and surface energy pa- 
rameters, we expect our results to be insensitive to variations 
in the parameters as long as they describe the same (R ^> t,, 
and strong anchoring). 

5 Conclusions 

In this paper, we studied the interaction of two spheri- 
cal colloidal particles with degenerate planar anchoring in 
a nematic media by numerically minimizing the Landau-de 
Gennes 2 - bulk and Fournier— surface energy using a finite 
element method. Our choice of parameters belong to the 
regime of large particles (in comparison to the nematic co- 



a 



b 



c 




Fig. 9 Director profile and topological defects on the surface of the 
colloidal particles when they approach each other along the nematic 
direction (9 = 0) at different inter-particle separations (a) 
d/R = 3.00, (b) d/R = 2.70, and (c) d/R = 2.10. The regions where 
the scalar order parameter drops below 0.65/, is highlighted, 
signaling the existence of a topological defect. The defect points on 
the surface and bulk are indicated by blue and blue colors, 
respectively. 
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herence length) and strong surface anchoring. We obtained 
the nematic-induced effective potential energy of the system 
for different inter-particle separations and orientations with re- 
spect to the bulk nematic director. 

By studying the free energy landscape of the system, we 
found that the system assumes its unique global minimum of 
energy when the particles are in close contact and are oriented 
at an angle 8 = 28° ± 2 with respect to the bulk nematic direc- 
tor. Our results are in a very close agreement with the exper- 
imental results in Ref. 4, ~ 30°. To the best of our knowl- 
edge, we have provided the first clear theoretical explanation 
of these experimental findings. 

Our results suggest that for large inter-particle separations, 
the quadrupolar structure of the defect-pair on each of the 
particles is essentially preserved, resulting in a monotonous 
attractive or repulsive inter-particle net force, depending on 
the orientation angle. However, for smaller orientation angles 
(9 > 15°) and at smaller inter-particle separations, the axial 
symmetry of the defect-pairs is continuously broken, resulting 
in the emergence of an attractive interaction due to the ten- 
dency of the system to reduce the volume of distorted fluid. 
This very unexpected attraction, in very short distances, has 
not been reported before and may be of interest to be explored 
experimentally too. 

The Finite Element method, used in this study, simply can 
be extended to more complicated geometries. Inter-particles 
interaction for nonspherical colloids with planar anchoring™ 
and also many-body interactions between the colloids in col- 




Fig. 11 Director profiles on the surface of the colloids for a center 
separation of d/R = 2. 10, and for four different orientations (a) 
6 = 20°, (b) 6 = 30°, (c) 6 = 40° and (d) 6 = 50°. The spheres are 
viewed along the bulk nematic director and from the side where the 
boojums are closer to each other. The center of each of the spheres 
is indicated by a blue dot. It is noticed that the axial symmetry of the 
defect-pair on each of the spheres is essentially preserved. 
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e=30 




Fig. 12 Defect displacement as a function of inter-particle 
distances for 9 = 0°,30°,60° and 90°. Displacements are scaled to 
the particles radius. It shows that defects do not moves a lot for large 
angles. The curves are only as eye guides. 



loidal aggregation s i are in direction of our future studies. 
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